function rslt = delta_r(Name, Data, Params, initial_rate)
% Name = 'R_noDriftRSCEVGARCH';
% load data;
S = Data;%num(1:R);
ret = S(2:end) - S(1:end-1);
R_t_2 = S(end);
h_t_1 = var(ret);
R_t_1 = initial_rate;
name = Name(3:end);

switch Name(1)
    case 'S'    %single factor
        switch(name)
            case 'lognormal'
                alpha_1 = Params(1);
                sigma = Params(2);
                mean = alpha_1.*R_t_1;
                volatility = sigma.*R_t_1;
            case 'dothan'
                mean = 0;
                sigma = Params(1);
                volatility = sigma.*R_t_1;
            case 'pureCEV'
                mean = 0;
                sigma = Params(1);
                rho = Params(2);
                volatility = sigma.*(R_t_1.^rho);
            case 'Vasicek'
                alpha_0 = Params(1);
                alpha_1 = Params(2);
                sigma = Params(3);
                mean = alpha_0 + alpha_1*R_t_1;
                volatility = sigma;
            case 'CIR'
                alpha_0 = Params(1);
                alpha_1 = Params(2);
                sigma = Params(3);
                mean = alpha_0 + alpha_1*R_t_1;
                volatility = sigma.*sqrt(R_t_1);
            case 'CKLS'
                alpha_0 = Params(1);
                alpha_1 = Params(2);
                sigma = Params(3);
                rho = Params(4);
                mean = alpha_0 + alpha_1*R_t_1;
                volatility = sigma.*(R_t_1.^rho);
            case 'nonlinear'
                alpha = Params(1);
                alpha_0 = Params(2);
                alpha_1 = Params(3);
                alpha_2 = Params(4);
                sigma = Params(5);
                rho = Params(6);
                mean = alpha./R_t_1 + alpha_0 + alpha_1.*R_t_1 + alpha_2.*R_t_1.^2;
                volatility = sigma.*(R_t_1.^rho);
            otherwise                       % RW by default
                alpha_0 = Params(1);
                sigma = Params(2);
                mean = alpha_0;
                volatility = sigma;
        end
        
    case 'G'    %GARCH
        beta_0 = Params(1);
        beta_1 = Params(2);
        beta_2 = Params(3);
        switch(name)
            case {'noDrift','noDriftCEV'}
                mean = 0;
                
            case {'linearDrift','linearDriftCEV'}
                alpha_0 = Params(4);
                alpha_1 = Params(5);
                mean = alpha_0 + alpha_1.*R_t_1;
                
            case {'nonlinearDrift','nonlinearDriftCEV'}
                alpha = Params(4);
                alpha_0 = Params(5);
                alpha_1 = Params(6);
                alpha_2 = Params(7);
                mean = alpha./R_t_1 + alpha_0 + alpha_1.*R_t_1 + alpha_2.*R_t_1.^2;
        end
        residual = R_t_2 - mean;
        h_t = beta_0 + beta_1.*residual.^2 + beta_2.*h_t_1;
        switch(name)
            case {'noDriftCEV','linearDriftCEV','nonlinearDriftCEV'}
                rho = Params(end);
                volatility = sqrt(h_t).*(R_t_1.^rho);
            otherwise
                volatility = sqrt(h_t);
        end
        
%     case 'R'    %regime switching
%         offset = 1;
%         
%         c1 = Params(offset);
%         d1 = Params(offset+1);
%         sigma_1 = 1;
%         
%         c2 = Params(offset+2);
%         d2 = Params(offset+3);
%         sigma_2 = Params(offset+4);
%         
%         p = 1./(1+ exp(-c1 - d1.*R_t_1));% p11 = P(s_t = 1 | s_{t-1} = 1)
%         q = 1./(1+ exp(-c2 - d2.*R_t_1));% p22 = P(s_t = 2 | s_{t-1} = 2)
%         
%         offset = offset + 5;
%         
%         switch(name)
%             
%             case {'noDriftRSCEV','noDriftRSGARCH','noDriftRSCEVGARCH'}
%                 mean1 = zeros(n,1);
%                 mean2 = zeros(n,1);
%                 
%             case {'linearDriftRSCEV','linearDriftRSGARCH','linearDriftRSCEVGARCH'}
%                 regime
%                 alpha_01 = Params(offset);
%                 alpha_11 = Params(offset+1);
%                 mean1 = alpha_01 + alpha_11.*R_t_1;
%                 
%                 regime 2
%                 alpha_02 = Params(offset+2);
%                 alpha_12 = Params(offset+3);
%                 offset = offset +4;
%                 mean2 = alpha_02 + alpha_12.*R_t_1;
%                 
%             case {'nonlinearDriftRSCEV','nonlinearDriftRSGARCH','nonlinearDriftRSCEVGARCH'}
%                 regime 1
%                 alpha_a1 = Params(offset);
%                 
%                 alpha_01 = Params(offset+1);
%                 alpha_11 = Params(offset+2);
%                 alpha_21 = Params(offset+3);
%                 mean1 = alpha_a1./R_t_1 + alpha_01 + alpha_11.*R_t_1 + alpha_21.*R_t_1.^2;
%                 
%                 regime 2;
%                 alpha_a2 = Params(offset+4);
%                 alpha_02 = Params(offset+5);
%                 alpha_12 = Params(offset+6);
%                 alpha_22 = Params(offset+7);
%                 
%                 mean2 = alpha_a2./R_t_1 + alpha_02 + alpha_12.*R_t_1 + alpha_22.*R_t_1.^2;
%                 
%                 offset = offset + 8;
%         end
%         residual_1 = ret - mean1;
%         residual_2 = ret - mean2;
%         xi = zeros(n,1);
%         xi_0 = 0.5;
%         h0 = 7.24814E-05;
%         switch(name)
%             case {'noDriftRSCEV','linearDriftRSCEV','nonlinearDriftRSCEV'}
%                 rho_1 = Params(offset);
%                 rho_2 = Params(offset+1);
%                 offset = offset + 2;
%                 volatility_1 = sigma_1.*(R_t_1.^rho_1);
%                 volatility_2 = sigma_2.*(R_t_1.^rho_2);
%                 eta_1 = normpdf(ret,mean1,volatility_1);
%                 eta_2 = normpdf(ret,mean2,volatility_2);
%                 if(cdf)
%                     fyCDF = zeros(n,1);
%                     eta_1CDF = normcdf(ret,mean1,volatility_1);
%                     eta_2CDF = normcdf(ret,mean2,volatility_2);
%                 end
%                         if( sum(eta_1 > 10000)||sum(eta_2>10000)|| all((eta_1==0)&(eta_2==0)) )
%                 if(all((eta_1==0)&(eta_2==0)))
%                     f = -10000;
%                     return;
%                 end
%                 fy = zeros(n,1);
%                 
%                 for i = 1: n
%                     
%                     if(i==1)
%                         xi_i_1 = xi_0;%xi{i-1}
%                     else
%                         xi_i_1 = xi(i-1);
%                     end
%                     fy(i) = p(i).*xi_i_1.*eta_1(i) + (1-p(i)).*xi_i_1.*eta_2(i) +...
%                         q(i).*(1-xi_i_1).*eta_2(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i);
%                     fyCDF(i) = p(i).*xi_i_1.*eta_1CDF(i) + (1-p(i)).*xi_i_1.*eta_2CDF(i) +...
%                         q(i).*(1-xi_i_1).*eta_2CDF(i) + (1-q(i)).*(1-xi_i_1).*eta_1CDF(i);
%                     xi(i) = (p(i).*xi_i_1.*eta_1(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i))./fy(i);
%                 end
%                 
%             case {'noDriftRSGARCH','linearDriftRSGARCH','nonlinearDriftRSGARCH'}
%                 beta_0 = Params(offset);
%                 beta_1 = Params(offset+1);
%                 beta_2 = Params(offset+2);
%                 offset = offset + 3;
%                 
%                 eta_1 = zeros(n,1);
%                 eta_2 = zeros(n,1);
%                 if(cdf)
%                     eta_1CDF = zeros(n,1);
%                     eta_2CDF = zeros(n,1);
%                     fyCDF = zeros(n,1);
%                 end
%                 e1 = zeros(n,1);
%                 e2 = zeros(n,1);
%                 h = zeros(n,1);
%                 volatility_1 = zeros(n,1);
%                 volatility_2 = zeros(n,1);
%                 fy = zeros(n,1);
%                 
%                 for i = 1: n
%                     if(i==1)
%                         xi_i_1 = xi_0;  %xi{i-1}
%                         h(i) = h0;
%                     else
%                         xi_i_1 = xi(i-1);
%                         e1(i) = (residual_1(i-1)./sigma_1);
%                         e2(i) = (residual_2(i-1)./sigma_2);
%                         h(i) = beta_0 + ...
%                             (e1(i).^2 .* xi_i_1 + e1(i).^2 .* (1-xi_i_1))*beta_1 +...
%                             beta_2 * h(i-1);
%                         
%                     end
%                     if(h(i)==0)
%                         f = -10000;
%                         return;
%                     end
%                     volatility_1(i) = sqrt(h(i))*sigma_1;
%                     volatility_2(i) = sqrt(h(i))*sigma_2;
%                     
%                     eta_1(i) = normpdf(ret(i),mean1(i),volatility_1(i));
%                     eta_2(i) = normpdf(ret(i),mean2(i),volatility_2(i));
%                     if(cdf)
%                         eta_1CDF(i) = normcdf(ret(i),mean1(i),volatility_1(i));
%                         eta_2CDF(i) = normcdf(ret(i),mean2(i),volatility_2(i));
%                         fyCDF(i) = p(i).*xi_i_1.*eta_1CDF(i) + (1-p(i)).*xi_i_1.*eta_2CDF(i) +...
%                             q(i).*(1-xi_i_1).*eta_2CDF(i) + (1-q(i)).*(1-xi_i_1).*eta_1CDF(i);
%                     end
%                     
%                                 if(eta_1(i) > 10000)|| eta_2(i)>10000 || ((eta_1(i)==0)&&(eta_2(i)==0))
%                     if(all((eta_1==0)&(eta_2==0)))
%                         f = -10000;
%                         return;
%                     end
%                     
%                     fy(i) = p(i).*xi_i_1.*eta_1(i) + (1-p(i)).*xi_i_1.*eta_2(i) +...
%                         q(i).*(1-xi_i_1).*eta_2(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i);
%                     xi(i) = (p(i).*xi_i_1.*eta_1(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i))./fy(i);
%                 end
%             case {'noDriftRSCEVGARCH','linearDriftRSCEVGARCH','nonlinearDriftRSCEVGARCH'}
%                 rho_1 = Params(offset);
%                 rho_2 = Params(offset+1);
%                 beta_0 = Params(offset+2);
%                 beta_1 = Params(offset+3);
%                 beta_2 = Params(offset+4);
%                 offset = offset + 5;
%                 if(cdf)
%                     eta_1CDF = zeros(n,1);
%                     eta_2CDF = zeros(n,1);
%                     fyCDF = zeros(n,1);
%                 end
%                 eta_1 = zeros(n,1);
%                 eta_2 = zeros(n,1);
%                 
%                 e1 = zeros(n,1);
%                 e2 = zeros(n,1);
%                 h = zeros(n,1);
%                 volatility_1 = zeros(n,1);
%                 volatility_2 = zeros(n,1);
%                 fy = zeros(n,1);
%                 
%                 for i = 1: n
%                     
%                     if(i==1)
%                         xi_i_1 = 0.5;%xi{i-1}
%                         h(i) = h0;
%                     else
%                         xi_i_1 = xi(i-1);
%                         e1(i) = (residual_1(i-1)./sigma_1);
%                         e2(i) = (residual_2(i-1)./sigma_2);
%                         if(h(i-1)==0)
%                             f = -10000;
%                             return
%                         end
%                         h(i) = beta_0 + ...
%                             (e1(i).^2 .* xi_i_1 + e1(i).^2 .* (1-xi_i_1))*beta_1 +...
%                             beta_2 * h(i-1);
%                         
%                     end
%                     if(h(i)==0)
%                         f = -10000;
%                         return
%                     end
%                     volatility_1(i) = (R_t_1(i).^rho_1).*sqrt(h(i))*sigma_1;
%                     volatility_2(i) = (R_t_1(i).^rho_2).*sqrt(h(i))*sigma_2;
%                     eta_1(i) = normpdf(ret(i),mean1(i),volatility_1(i));
%                     eta_2(i) = normpdf(ret(i),mean2(i),volatility_2(i));
%                     
%                     if(cdf)
%                         eta_1CDF(i) = normcdf(ret(i),mean1(i),volatility_1(i));
%                         eta_2CDF(i) = normcdf(ret(i),mean2(i),volatility_2(i));
%                         fyCDF(i) = p(i).*xi_i_1.*eta_1CDF(i) + (1-p(i)).*xi_i_1.*eta_2CDF(i) +...
%                             q(i).*(1-xi_i_1).*eta_2CDF(i) + (1-q(i)).*(1-xi_i_1).*eta_1CDF(i);
%                     end
%                     
%                                 if(eta_1(i) > 10000 ||eta_2(i)>10000 || (eta_1(i)==0 &&eta_2(i)==0))
%                     if(all((eta_1==0)&(eta_2==0)))
%                         f = -10000;
%                         return;
%                     end
%                     
%                     fy(i) = p(i).*xi_i_1.*eta_1(i) + (1-p(i)).*xi_i_1.*eta_2(i) +...
%                         q(i).*(1-xi_i_1).*eta_2(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i);
%                     
%                     xi(i) = (p(i).*xi_i_1.*eta_1(i) + (1-q(i)).*(1-xi_i_1).*eta_1(i))./fy(i);
%                     
%                 end
%         end
%         
%     case 'J'    %jump diffusion
%         c = Params(1);
%         d = Params(2);
%         q = 1./(1+exp(-c-d.*R_t_1));
%         offset = 3;
%         meanX = 0;
%         switch(name)
%             case {'noDriftJDCEV','noDriftJDGARCH','noDriftJDCEVGARCH'}
%                 meanX = 0;
%                 
%             case {'linearDriftJDCEV','linearDriftJDGARCH','linearDriftJDCEVGARCH'}
%                 alpha_0 = Params(offset);
%                 alpha_1 = Params(offset+1);
%                 meanX = alpha_0 + alpha_1.*R_t_1;
%                 offset = offset + 2;
%                 
%             case {'nonlinearDriftJDCEV','nonlinearDriftJDGARCH','nonlinearDriftJDCEVGARCH'}
%                 alpha = Params(offset);
%                 alpha_0 = Params(offset+1);
%                 alpha_1 = Params(offset+2);
%                 alpha_2 = Params(offset+3);
%                 offset = offset + 4;
%                 meanX = alpha./R_t_1 + alpha_0 + alpha_1.*R_t_1 + alpha_2.*R_t_1.^2;
%         end
%         sigma = Params(offset);
%         muy = Params(offset+1);
%         gamma = Params(offset+2);
%         offset = offset + 3;
%         
%         n = length(ret);
%         compound_meanX = meanX + q.*muy;
%         
%         
%         
%         switch(name)
%             case {'noDriftJDCEV','linearDriftJDCEV','nonlinearDriftJDCEV'}
%                 rho = Params(offset);
%                 volatility = sigma.*(R_t_1.^rho);
%             case {'noDriftJDGARCH','linearDriftJDGARCH','nonlinearDriftJDGARCH'}
%                 beta_0 = Params(offset);
%                 beta_1 = Params(offset+1);
%                 beta_2 = Params(offset+2);
%                         offset = offset + 3;
%                 if(beta_0 < 0 || beta_1 < 0 || beta_2 < 0)
%                     f = -1000000;
%                     return;
%                 end
%                 
%                 h = zeros(n,1);  % 2 -> N
%                 residual = ret - compound_meanX;
%                 h(1) = var(ret); %unconditonal variance
%                 for i = 2: n
%                     h(i) = beta_0 + beta_1.*residual(i-1).^2 + beta_2.*h(i-1);
%                 end
%                 volatility = sigma.*sqrt(h);
%                 
%             case {'noDriftJDCEVGARCH','linearDriftJDCEVGARCH','nonlinearDriftJDCEVGARCH'}
%                 beta_0 = Params(offset);
%                 beta_1 = Params(offset+1);
%                 beta_2 = Params(offset+2);
%                 if(beta_0 < 0 || beta_1 < 0 || beta_2 < 0)
%                     f = -1000000;
%                     return;
%                 end
%                 rho = Params(offset + 3);
%                         offset = offset + 4;
%                 h = zeros(n,1);  % 2 -> N
%                 residual = ret - compound_meanX;
%                 h(1) = var(ret); %unconditonal variance
%                 for i = 2: n
%                     h(i) = beta_0 + beta_1.*residual(i-1).^2 + beta_2.*h(i-1);
%                 end
%                 volatility = sigma.*sqrt(h).*(R_t_1.^rho);
%         end
end

rslt = mean + volatility*normrnd(0,1);